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ABSTRACT 

The  neutron  transport  equation  is  differenced  with 
respect  to  all  but  the  spatial  variables  and  the  solutions 
of  the  resulting  systems  of  ordinary  differential  equations 
are  studied.   For  the  case  of  elastic  scattering  the  approxi- 
mate solution  is  obtained  and  its  convergence  to  the  exact 
solution  is  proven.   If  very  accurate  approximations  are 
required  the  solution  may  be  evaluated  by  computing  machines. 
It  is  seen  that  the  present  procedure  has  many  advantages 
over  direct  niimerical  integration  procedures. 

For  elastic  scattering  the  procedure  is  just  the  well- 
known  Wick-Chandrasekhar  method.   The  solution  for  homogeneous 
slabs  of  finite  thickness  is  reduced  to  the  inversion  of  at 
most  three  matrices  of  small  order.   The  required  calculations 
are  easily  done  on  desk  computers,   A  new  relation  between 
escape  probability  and  capture  fraction  is  also  obtained. 
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APPROXIMATE   SOLUTIONS   OP   STEADY-STATE 
NEUTRON  TRANSPORT    PROBLEMS  FOR   SLABS 

1«   Introduction. 

In  this  paper  we  study  some  methods  for  the  approximate 
solution  of  various  steady  state  neutron  transport  problems. 
The  methods  consist  in  replacing  some,  but  not  all,  of  the 
continuous  independent  variables  by  a  discrete  set  of  mesh 
points  and  then  differencing   the  transport  equation  with 
respect  to  these  discrete  variables.   An  approximating  system 
of  functional  equations  is  obtained  In  this  manner.   In  the 
problems  treated  in  this  paper  the  approximating  equations 
are  systems  of  ordinary  differential  equations.   Such 
procedures  have  been  used  in  [2,9]  for  unsteady  problems  in 
which  case  the  approximating  equations  form  hyperbolic 
systems  of  partial  differential  equations.   The  approximating 
procedures  of  this  paper  are  closely  related  to  methods  which 
have  been  proposed  for  the  direct  numerical  solution  of 
transport  problems  [l,2,3»l4-,5]  • 


The  term  "differencing"  is  used  in  a  generalized  sense 
which  may  include  integrating  or  averaging  the  equations 
over  a  mesh  interval  as  well  as  the  use  of  numerical 
integration  techniques. 
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In  Part  I  (Sections  2-5)  we  consider  the  elastic  scatter- 
ing transport  equation  for  an  infinite  plane  slab  of  finite 
thickness  and  arbitrary  composition.   The  convergence  of  the 
approximate  solution,  which  is  obtained  explicitly,  to  the 
exact  solution  in  the  appropriate  limit  is  demonstrated© 
It  is  also  shown  that  the  approximate  solution  is  stable 
with  respect  to  the  growth  of  small  errors  introduced  by 
approximation  or  inaccurate  knowledge  of  the  inhomogeneous 
source  terms.   The  method  and  solution  are  easily  modified 
to  solve  the  spherically  symmetric  case. 

In  Part  II  (Sections  6-10)  we  consider  the  isotropic 
scattering  mono-energetic  transport  equation  for  homogeneous 
slabs.   The  method  for  such  problems  is  precisely  the 
Wick-Ghandrasekhar  procedure  [10,11].   While  this  method  is 
discussed  at  some  length  in  the  literature  it  is  believed 
that  the  present  discussion  simplifies  its  application  to  a 
variety  of  problems  concerned  with  slabs  of  finite  thickness. 
The  problems  discussed  are  the  determination  of  criticallty, 
capture  fraction  and  escape  probability.   The  convergence  of 
the  method  is  not  proved  but  a  new  result  relating  escape 
probability  and  capture  fraction  is  obtained. 
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PART  I 

Elastic  Scattering 

2»     Formulation  of  Elastlc-scatterlng    Problems. 

The   steady  state  neutron  transport   equation  In  plane 
geometry   is 

(2.0a)  T[i(x,ti,u)]   =p  g^  +  (5;j,(x,u)|  i(x,n,u)-S[5||x,u]   = 

Here  J     2(x,(j,,u)      is   the  flux   of  neutrons   at  position     x, 

■with  lethargy     u,      whose  velocity  vector  makes   an  angle 

cos"     \i     with  the   positive     x-axis  |      cJ^(x,u)      is    the   total 

macroscopic    cross-section  at      x     for  neutrons  with  lethargy 

u|      ^(x,ti.,u)      is   an   inhomogeneous   source   of  neutrons^     S      is 

the  scattering   term.      For   elastic   scattering    this    term  can  be 

written  as 

u 

S[J|ii,u]   =      J      g(xfu,u')   H[||p.,u,u«]   du«      , 
u-q 
(2o0b) 


n 


H[I|l^,u,u«  ]   =    J  i(x,m(ti,u-u«,0),u' )   de      . 

o 

E'or  a  scattering  element  with  atomic  weight  A  the  quanti- 
ties to  be  used  in  (2,0b)  are  defined  ass 

rq=^n(f±i)',   gUju,u.)  ^c^(x,uM  ^^e-'-^   , 

(2.1)  <  m(ti,y,e)  =  ^a{j)   +   [l-ti^l^^^  [l-a^ij)]^^^   cos  6   , 

(^  a(y)  =  i  [(A+De^/^  „  (^_3^)  ^-y/2j   ^ 
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Here      G^  (x,u)      is    the  macroscopic   scattering  cross-section 
at     X     for  neutrons   with  lethargy     u.      The   derivation  and 
physical    interpretation   of   the   quantities    in    (2,1)    are 
thoroughly  discussed   in   the   literature    [ 5,6,7 ]•      If  more 
than   one   scattering'  element   is   present   additional   scatter- 
ing  terras    of   the  form    (2.0b)   must  be   included   in    (2o0a), 
In  spherically  symmetric   geometry  the    scattering    terra 
retains    the  above  form  while   the  differential   operator   in 
(2,0a)    becomes 

Such  problems  may  be  treated  by  slightly  modifying  the 
present  method,  say  in  a  manner  clearly  shown  in  [1]  or  [2], 

The  flux  of  neutrons  entering  the  slab,   -a  <  x  <  a, 
from  the  outside  must  be  specified.   If  we  assume  no  sources 
present  in  the  vacuum  external  to  this  slab  there  can  be  no 
neutrons  entering  it  and  we  have  the  boundary  conditions 

(2.2)   J(-a,tt,u)  =0,    0  <  ti  <  1  J   i(a,ii,u)  =0,    ~1  <   \i  <   0, 

A  specified  non-zero  flux  incident  on  the  slab  can  be 
included  in  the  present  formulation  either  by  modifying 
(2,2)  to  include  Inhomogeneous  terms  or  by  including 
spatial  Dirac  5-f unctions  in  the  inhomogeneous  source 
J   (x,ii,u)   in  (2.0a). 

Since  only   elastic  scattering  is  considered  there 
can  be  no  neutrons  present  with  energy  greater  than  that 
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of  the  maxlmtim  energy  source  neutrons,   (If  fission  were 
allowed  the  maximum  energy  would  be'  that  of  the  fission 
spectrum.)   If  this  energy  is  E  ;y  and  the  lethargy 
variable  corresponding  to  any  energy  E  is  defined  by 

E 


u  =  in 


max 
E    » 


then     u  >   0     for  all  neutrons   In  the   present  problem.     That 
is     i(x,^L,u)   =  0     for     u  <   0     and  by    (2,0b),      S[5jp.,0]   =  0. 
Thus   at   the  maximum  energy    (2,0a)   simplifies   and  can  be 
solved  subject   to   the   bovindary   conditions    (2,2),      The 
solution  is 

X  X 

\^li(^.t^»0)    exp[^  J   0-^(4', 0)d^']dC,    0<t.<l| 
(2.3)   i(x,ti,0)    =< 


-^f    J(C,M.,0)    exp[f^  J  cr^(C',0)d?']dC,    -l<ti<0. 


If  fission   is   present   or    the   scattering   is    isotropic 
SEjjM-fO]   ^  0     and   the   flux   of  maximum  energy  neutrons 
cannot  be   simply   obtained  as   above. 
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3»  Reduction  to  Differential-Difference  Equations. 

The  continuous  variables   (lJ.,u)   are  to  be  replaced  by 
a  mesh  of  discrete  points   (|x.,u  )   at  each  of  which  the  flux, 
5(x,|j-.,u  ),   is  to  be  approximated  by  a  function  ^(x). 
These  approximating  functions  are  determined  as  a  solution 
of  a  system  of  linear  differential  equations  which  in  some 
sense  approximate  the  transport  equation  (2,0) »    There  is 
no  unique  procedure  by  which  to  derive  systems  of  approximat- 
ing equations  and  a  priori  it  cannot  be  determined  which  of  a 
class  of  "equivalent  systems"  will  yield  the  best  approximations 
to  J.   (Two  systems  are  equivalent  if  they  approximate  (2,0)  to 
the  same  order  of  magnitude  in  the  mesh  spacing »)   Of  course  any 
approximating  system  should  satisfy  the  conditions  of 
consistency  (approximating  equations  — *•  transport  equation  as 

Aix,  Au  — s»  0),   convergence   (i^j(x)  -~**  l^^yl^jj^n^   ^^ 
AM-,  z\u  — >   0)   and  stability  ("small"  changes  in  Inhomogeneous 
terms  produce  "small"  changes  in  the  solution).   These  condi- 
tions are  defined  more  precisely  in  Section  5  aJ^d  are  shown 
to  hold  for  the  approximating  system  derived  in  the  present 
section. 

Since,  in  actual  computations,  the  limit  as  aix,  au  — >  0 
cannot  be  attained  additional  conditions  are  frequently 
imposed  on  the  approximating  equations  to  insure  that  their 
solutions  have  some  required  physical  properties o   For 
instance,  since  J  is  a  neutron  flux,  it  is  reasonable  to 
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require  that  )^(x)  >  0  or,  as  Carlson  [l]  does,  to  require 
that  the  total  number  of  nuetrons  be  conserved.   In  the 
present  paper  only  the  positivity  condition  will  be 
consideredo 

The  approximating  equations  will  be  derived  in  two 
steps,  first  differencing  with  respect  to  the  angular 
variable,  \i,      and  second  differencing  with  respect  to 
the  lethargy  variable,   Ue   This  procedure  will  exhibit 
two  different  (and  not  equivalent)  approximating  systems  of 
functional  equations.   It  is  then  easy  to  modify  the  procedure 
and  a  clear  physical  interpretation  of  the  equations  in  each 
stage  is  made  obvious o   The  mesh  spacing  is  taken  to  be 
uniform  only  for  simplicity  of  presentation. 
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3A.   Angular  Differencing. 

The  ^.-Interval,   -1  <  t^  <  1,   is  divided  Into   (2J-1) 
sublntervals  of  length  A  p.  by  the  mesh  points 

(3oO)    ^.  =  -^  (|j|  -  i)  Aix,   l<|jl<J,   ^M-  =  pfe  • 

This  mesh  has  been  chosen  to  avoid   ^i  =  0  as  a  mesh  point. 
Such  a  point  could  be  included  provided  the  equations 
corresponding  to  (3o2)  in  the  new  mesh  are  "centered"  at 

^^1+1/2  ^^  ^^  ^°^®  ^^  ^^^    ^^^  '■^^'  '^^^   flux  and  scattering 
term  at  each  angular  mesh  point  are  to  be  approximated  by 
the  quantities 

(3.1)  i(x,ttj,u)  ^  ^j(x,u)  ,    S[ijtij,u]  ^  S^ijU^iu]    . 

Replacing   corresponding   terms    In    (2,0a)    by   the   above 
approximations   and   setting      [i  =  |J. .      we   obtain  for    each      j: 

J 

(3.2)  T:^^[0.{x,n)]   =  jttj  ^  +  o^(x,u)|^^(x,u)-Sj.[jZ(j^ju]   =  ;/j  (x,u), 

1  <    Ul    <  J   . 

Here  we   have   used   the   obvious  notation       7/.  (x,u)    =  3j(x,ii  .,u) , 
By  specifying   the   dependence   of   the   approximate   scattering 
term,      S>.[0-^,u],      on   the   approximate  fluxes,      0^{x,u) ,       (3<>2) 
becomes    a  system   of  functional    equations   in   the  imknowns 
^.(x,u) , 

The  form   of     S.      Is   determined  by  using    a  nvunerlcal 
Integration   technique   in    (2. Ob),      A  straightforward  procedure 

-   11  - 


is  to  approximate  the  integrand  in  H  by  a  piecewise  linear 
function: 


(3.3) 


when 


m(ix.,y,e)-(x, 

^-^  <  y  <  ^  • 


Integrating    this   approximate    integrand  with  respect   to     Q 

over   the   interval      0  <   Q  <  %     we   obtain,    provided      0  <  u-u'    <  q, 

J 
(3.i|.)      H^.[jZJjj.;u,u' ]   s     211     ^ji^(u-u')^lj.(x,u')    ^  H[ijM.j,u,u»]    , 

|k|=l 

The  weights,        '^ik'      '"^^  ^®  derived  by  carrying   out   the 

indicated    integrations*     We   do  not   include   them  here  but 

rather   list  their  most   important  properties? 

J 
(3.5)      U)jk(y)   >  0»       m  ^jk^y^   "  ^»      'ijk^^^   "  ^^jk»    ^^SJ<!=i)o 

|k|=l 

(The  last  of  these  properties  is  crucial  in  obtaining  the 

simple  explicit  solution  of  Section  i|.o)   If  a  more  accurate 

interpolation  is  used  in  place  of  (3o3)  different  weights 

are  obtained  in  (3.i;)»   However,  these  weights  still  satisfy 

the    last    two  conditions    of    (3«5)    and   can  be  made  non-negative 

by  requiring  that  the  interpolation  formula  be  essentially  of 

the  closed  type  (see  [8]),  We  assume  then  an  approximation 
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of  the  form  O.I^)  with  weights  satisfying  (3»5)  and  an 
accuracy  not  less  than  woiild  be  obtained  by  using  linear 
interpclaticn  as  in  {3«3)»   -he  aprroxinate  scattering 
term  is  now  defined  by  using  {3»L)   in  (2.0b)  to  get 

S..[j2j^;u]  =  J  g(xju,u'  )  H.[^,  ;u,u']  du«   , 


J 


(3.6)  J    u 

=  YZ       I  g(xju,u»)  :^  (u-u'  )  ^^(x,u' )  du'. 

|kl=l  '^'^ 

Equations  (3«2)  and  ,3»c)  fcm  a  srsterr.  of   2J 
integro -differential  equations,  the  first  approximating 
system  J  for  the  determination  of  the  ^,.  (x,u),   the  first 
approximating  fl^jxes.   The  boundary  conditions  to  be 
imposed  are,  from  (2,2.), 

(3.7)  )^.(-a,u)  =  C,   -  <  :"  <  -'»    ^j(a,u)  =  C,   -J  <  j  < -L 

As  in  Section  2  we  may  solve  for  ^.(x,0)   since   3.[^,  ,0]  =  0. 
1r.e    solution  is  exactly  of  the  form  (2.3)  and  we  find  that 

(3.6)  /5.(x,C)  =  i(x,^..,0)  . 

J  J 


3B>   Lethargy  Differencing. 

The  lethargy  Interval,   0  <  u  <  U,   Is  divided  into 
N  Intervals  of  length   Au  by  the  mesh  points 

q    1 
(3«9)     u^  =  nAu,   0  <  n  <  N,   Au  =  ^  (=  ^  if   q  =  oo  )o 

With  this  choice  of  spacing  the  lethargy  integrations  in 
(3»6)  are  over  M  Intervals  (provided  u  >  q).   If  more 
than  one  scattering  element  is  present  the  integrations 
in  the  additional  scattering  terms  may  be  over  a  non-integral 
number  of  u-intervals,,   However  an  interpolation  procedure, 
as  Indicated  in  [^],  could  be  used  to  approximate  the 
Integrals  over  fractional  lethargy  Intervals.   In  practical 
computations  the  mesh  (3.9)  is  likely  to  be  chosen  as 
non-uniforra  to  eliminate  the  need  for  many  lethargy  points 
if  u/q  »  1  for  some  scatterer. 

We  now  approximate  the  first  approximating  fluxes  and 
scattering  terms  at  each  lethargy  mesh  point  by 

(3.10)  0^{x,u^)    ^iZ(^(x)  ,      Sj.[j2fj^ju^]  ^S^LiZJ^]  . 

Replacing  corresponding  teiTns  in  (3*2)  by  the  above 
approximations  we  obtain,  with  u  =  u  , 

(3.11)  T^;;  [iZl^(x)]  ^[t^j  ^-^  cr^^(x)jjzJ^(x)-s5[iZ(2]  =  ^^(x). 
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Here  we  have  used  oi,  (x)  s  o"m(x,u  )   and  y^^(x)  =5;<^(x,n.,u  ), 

The  explicit  form  of   S.  is  determined  by  using  numerical 

integration  in  (3«6).   If  the  trapezoidal  rule  is  used  we 
obtain 

n      J      nm         / 

where; 


nm 


,nm  /  „  \    .n-ra  ,   ^  _  ^  _  <iu 


-rLik(x)  =  V   (^>^k   »   ^o  =  ^n  =  T^'  ^m  =  ^^'  "^°'^J 


n-m 


(3.13)   «*Jjij.  5  £Ojj^([n-in]  Au)  } 


.2 


nm/  \    /       \  _  ^Bi/  \  (A+1)   -(n-m)AU 
g   (x)  =  g(x,u^,u^)  =cr^(x)  ^^  e 

If  a  more  accurate  integration  formula  is  used  (3.12)  and 
(3.13)  are  modified  only  by  a  change  in  the  definitions 
of  the   e  .   However,  by  again  requiring  a  formula  of  closed 
type  the  new  weights  are  still  positive.   For  lethargy  points 
u  <  q,   n-M  is  negative.   However,  since  there  are  no 
neutrons  with  negative  lethargy,  we  set  i^T,(x)  =  0  for  m  <  0 
and  in  these  cases  the  siim  (3.12)  starts  at  m  =  0.   If 
hydrogen  is  the  scatterer  q  =  oo   and  the  sum  always  starts 
from  m  =  0. 

The  system  of  differential  equations  (3»11-«12)  is  the 
second  approximating  system  for  the  determination  of  the  2NJ 
approximating  flxixes,   jZJ.(x)«   The  boundary  conditions  are, 
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by   (2.2)   or   (3.7), 

(3.11+)     fi^A-a)    =0,      1  <    j  5  JS     )Z^^(a)   =  0,    -J  <   j   <   -l;   0<n<N. 

Again  we  may    solve  for  the  ina:xiinuin   energy  flux,      ^.(x), 
since      S°[jZfJ]    =0.     We   find  that 

(3.15)  ^j(x)   =  jZJj(x,0)   =i(x,iij,0)    . 
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k»      Solution  of  Differential -Difference  Equations « 

The  system  for  the  determination  of  the  second  approxi- 
mating fluxes,  j^^(x),   can  be  solved  exactly  In  simple  closed 
form.   The  key  to  this  solution  Is  the  observation,  from  (3.5) 
and  (3.13),  that 

nn         -n-  Tin 
(i^.O)         -TLjk^^^  =  ^^  I  S^^J  ^jk  • 

Then   (3ol2)   can  be   written  as 

Sjt^k^    =     Au  I  g"^(x)   iZl^(x)    +  S^CjZf^]    , 

(i+.l)  r  -) 

^     «,  n-1       I       J  nm  ^         ( 

i  ■"="-«  l^|k|=i       '"^  j 

Using  this  result  In  (3»11)  yields 

ik.2)  [.J  ^  *   a-(.)  j  ^^(.)  =  l-l/^l   *  J-^.)   .r  ^  J^  ^  '' 

Here  a  modified  cross-section  has  been  Introduced  by  the 
definition 

a  (x)  =  (S^  (x)  -  Au  ^  g   (x)  , 

(l;.3) 

<5^  (x) 

Prom  (ij-.l)  It  Is  seen  that   S^  Is  Independent  of  j^^(x). 
Hence  for  each   (n,j)   (Ij.o2)  Is  a  linear  first  order 
Inhomogeneous  differential  equation  for  0.{x),      Bovmdary 
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conditions  are  given,  for  j  >  1  at  x  =  -a  and  for 
j  <  -1  at  X  =  a,   in  (3«lU)»   The  solutions  are  found 
to  be 

X  X 


ik.k)   0j(x)  = 


;^   f[S^[iZ^(5)]+i^(5)]  exp  [^  J"  a^(C')d4« 
^i     -a   ^  ^       ^         ^'j  C 

a  C 

^J  [sM(€)]+i'J(C)]  exp  [-i.  /  a^(C«)dC']d4, 


-J  <  j  <  -1  . 

starting  from  the  maximum  energy  flux,  i^j(x)   (which  is 
given  in  (3«l5))>  we  may  obtain  the  solution  for  all  n  by 
repeated  application  of  (i^-.l)  and  (i;,I|.),   If  the  source  and 
medium  are  piecewise  constant  all  integrations  can  be 
performed  explicitly.   The  solution  will  then  be  composed 
of  linear  combinations  of  exponentials  of  the  form! 

exp  (-^^''Vltijl). 

The  solution  (i]..I|.)  yields  only  positive  fluxes  (since 
all  cross-sections  and  the  source  are  non-negative).   However, 
the  flux  coming  directly  from  the  source  at  any  point  may  be 
exponentially  growing  if  a^(x)   is  negative  in  a  neighborhood 
of  that  point o   This  will  not  be  the  case  if,  as  seen  from 
(l4-«3)»  AU  is  sufficiently  small.   On  a  sequence  of  meshes 
in  which   au  — >  0   the   a  (x)   eventually  become  and  stay 
positive.   In  actual  computations,  however,  a  fixed  mesh  is 
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2 
used  and  the  above  considerations  suggest  that  we  require 


(i4-«5)  Au  <  max 

-a<x<a 


^05,^(X)      Q^ 


^J^i^)      (A+l)2. 


0<n<N 

With  such  a  choice  for   Au  the  exponential  decay  e'^'^'  '^  , 
rather  than  e      ''^  ,  may  be  considered  as  a  correction 
for  the  fact  that  some  neutrons  are  "scattered"  straight 
ahead  in  elastic  scattering  and  thus  the  total  cross-section, 
<yLf      is  too  large.   The  best  choice  for   Au  might  be 
determined  by  examining  some  experimental  data  which  measures 
this  foreward  scattering  effect*, 

In  the  case  of  spherically  symmetric  geometry  explicit 
solutions  can  also  be  obtained.   The  equations  corresponding 
to  (ij.«2)  will  be  coupled  pairwise  in  the  angular  subscript 
(that  is  0.      and  0 ^_-i      both  appear)  in  all  but  one  equation, 
say  that  for   j  =  J.   Starting  with  this  equation  the  system 
may  be  solved  recursively.   The  solution  is  of  the  form  (U.U) 
with  the  S.   replaced  by  successively  more  complicated 
expressions  as   j  decreases. 

If  many  angular  or  lethargy  mesh  points  are  used  the 
solutions  (i;»l4.)  may  become  too  complicated  for  direct  analysis. 

In  such  cases  a  numerical  evaluation  and  tabulation  of  the 

_  ^  ^  ________ 

This  condition  was  found  to  be  sufficient  in  [5]  for  the 
solvability  of  the  difference  equations  considered  there* 
With  an  additional  restriction  on  the  size  of  the  spatial 
mesh  it  also  guaranteed  positive  numerical  fluxes  and 
convergence  of  the  difference  solution. 
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approximate  solution  is  easily  accomplished  with  the  aid 
of  a  computing  machine,,   These  calculations  are  In  some 
sense  equivalent  to  direct  numerical  integrations  of  the 
transport  equationo   However,  the^ problem  of  evaluating 
expressions  of  the  form  occurring  in  (i4.oi|)  is  straight- 
forward and  thus  should  have  many  advantages  over  previously 
proposed  direct  numerical  procedures  [1,5] o 
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5«      Consistency,    Convergence   and  Stability. 

The   approximating   system    (3»11)   Is   consistent  with 
the    transport    equation    (2«0a)   if 


(5.0)        11m 

AU->0 


T[H/(x,ti,u)]    -  T^|^[^lf(x,tij.,u^)] 


I5IJI5J 


=   0,' 


\i=\j. .  0<  n  <N , 


u=u 


n 


for  all   sufficiently  smooth  functions      \|;(x,p.,u)«     To   demonstrate 
this   property  we  note   that,  by    (2«0),    (3o2),    (3«6),    (3ell)    and 
(3.12), 


,^Uf 


,AUr 


T[^l;]-T   ,'^[11;]    =    (T[>1/] -T^    [  ij,]  )    +    (T^,  U]    -  T^'^l^]) 


(S[\}f]    -   Sj[\l;])    +    (Sj[\l/]    -    S^[\lr]) 


u 


n 


It 


=     J       du'g(x,u^,u' )]  i  \|;(x,in(tx.,u^-u',9),u»)de  - 


u^-q 


^,.15.(%-U'  )     \lf(x,H^,U'   ) 


k  =1 


u 


n 


lkl=l/   %-^ 


J        du»g(x,u^,u'  )  ^,-ij(u^-u'  )\l;(x,iij^,u'  )- 


n  nm 

%  -n-ik(^>  >i/(x,ti^,u^) 


m=n-M 


Here  we  have  abbreviated  the  previous  notation  in  an  obvious 
manner.   Each  of  the  curly  brackets  above  contains  the  differ- 
ence between  some  Integral  and  an  approximation  to  that 


-  21  - 


integral*   These  differences  can  be  estimated  by  standard 
formulae •   If  the  simplest  procedures  described  in  Section  3 
are  used  in  these  approximate  integrations  it  is  easily 
shown  that 


(5.1)   |T[^]-T^j;[^]|  <  (aix)^  ^ 


+  (au)'=^  ^ 


Here  the  maxima  of  the  absolute  values  are  to  be  used.   If 

more  accurate  integration  formulae  had  been  used  the  right 

hand  side  of  (5ol)  would  contain  higher  powers  of   A  p,   and 

A  u»   It  is  clear  that  (5oO)  is  satisfied  provided  (since 

J  =  —  "*■  p)   the  passage  to  the  limit  is  such  that  au  /a\i   — >  0, 

Thus  the  equations  are  consistent  if  the  mesh  spacings  are 

related  by 


(5o2) 


A  u  <  cons 


t  (A.)'^/^-^^' 


for  any   e  >  0 


This  is  not  a  severe  restriction  and  could  be  relaxed  if  a 

3 
more  accurate  lethargy  integration  is  used  •   It  will  be 

shown  below  that  (5«2)  is  also  a  sufficient  condition  for 

convergence  and  stability. 

The  approximate  solution  will  converge  to  the  exact 

solution  provided 

If  the  error  in  the  lethargy  integration  is  of  order 

(/iu)^  condition  (5»2)  is  replaced  by;   au  <  const (a^,)^  /^  ^' 
for  any   e  >  0* 
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lim      li(x,ix.,u    )   -  0'\U)\    =  0    , 


1  <   hi  <  J 


0   <      n     <   N 


zsu->0 


To   prove   this  we   obtain   and   solve   a   system  of    equations 
for   the  error, 

n  n 


(5.3) 


e.(x)  =  i(x,iij,u^)   -  ^j(x) 


and  show  that  it  can  be  made  arbitrarily  small  as  the  mesh 

is  refined.   Setting  ^l  =  ix.   and   u  =  u   in  (2«0a)  and 

subtracting  (3«11)  from  it  yields,  using  (5.3), 

fl  1  1  j  I  <  J 
iS.k)       Oe^(x)]  =  S[J;,j,uJ  -  S-[i]  .  r^(x),{^  ^  ^ 


<  N  . 


Exactly  as  in  (5.1)  we  see  that 

2  1       '^'i 


{^.^)      lr^(x)|  <  A^^  ^ 


c)^ 


..u^^ 


-^(g^i) 


au 


=  r-. 


At  the  boundaries  the  error  satisfies,  from  (2«2)  and  (3<»7), 


(5.6)   ej(-a)  =0,   1  1  j  1  J  ;   e^(a)  =  0  ,   -J  5  j  <  -1  J 

0  <  n  <  N  • 

The  error  is  thus  determined  as  the  solution  of  the  system 
(5«l4-)  subject  to  the  boundary  conditions  (5.6 )»   This  system 
is  formally  identical  with  that,  (3«11)  and  (3.114-),  satisfied 
by  the  fluxes  jli^ix)      if  in  the  latter  we  replace  ^  .(x)  by 
tj^Cx).   Thus  by  the  procedure  of  Section  l\.   we  obtain  for  the 

J 


error 
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(5o7)  e^(x)   = 


X  X 

^      [[S^Ce^COl  +  tr^CC)]    exp[^  J  a^(?'  )dV  ]d^, 

a  ? 

^-   r  [S^[e^(4)]  +  T^(4)]exp[-^Ja^(?«)dC']d4, 

1  .      y  J       *^  J  U.  .       X 


m-i  X 


\  -J  <  j  <  -1  . 

To  estimate  the  magnitude  of  this  solution  we  Introduce  the 
quantities 

(5«8)  E  =   max  |e??(x)l  »  o"  =   max  cTg^Cx)  ,   as   mln  a^x)o 
-a<x<a  -a<x<a  -a<x<a 


-J<k<J 


-a<x<a 
0<n<N 


0<n<N 


nm 


Prom  (l|ol)  we  have,  since  .TL^'i^Ci)  >  0 

n-1 

\ r^ 


J      nm 
J  m=n-M   u  .  ^ 


.ni/ 


e-(4) 


(5.9) 


n-1     „^    )   J   n»m 
<  n  e^g^(C)<  T~  00 
m=n-M         ,  ,^.  3_ 


jk 


.m/ 


e»(«) 


n-1 


<  Tl 


m=n-M 


e  g'^(£)  E   , 
m^  ^^'     m  ' 


/ «  .T  \2  .   n-1 

^      ^        ^     (A+1)  ^-nAU  r; —   m^u  .p   _  ■^ 

<  Au  O^  ■  I. .    e      > e     £.  s=  b 

s   i^A  ^        m    n 


Here  we  have  used  the  definitions  (5o8),  (3ol3)  and  the 
properties  (3»5)«   Taking  the  absolute  value  In  (5o7)  and 
using  (5.9),  (5o8)  and  (5c5)  yields 
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(5.10)  |e^(x)|  <  [\-^^]  j'      a 

-~-J  exp[-=^(^-x)]d?,  -J  <  j  1  -1  . 

The  Integrals  may  be  evaluated  and  give,  respectively 

T     -a(a+x)/ti.  -,     -a(a-x)/lii.| 

i  (1-e        ^)    and    ^  (1-e  J  )   . 

We  note  that  these  quantities  are  positive  in  -a  <  x  <  a 
regardless  of  the  sign  of  a.  Then  taking  the  maximuin  in 
(5.10)  with  respect  to   x   and   j   yields 

if   a  >  0  , 


a 
The  above  inequality  can  be  written  as,  using  (5«9), 

(5.12)  e^^^E^  <      ^uG  Zle^^^E^  -.e^-^BT  , 

in=0 

where 

(5.13)  G  ^cr^i^^B      . 

Starting  with  n  =  1   and  applying  (5.12)  recursively 
yields  the  bound 

(5.11+)     E   <  Au[C  e  ""  (1+AuC)^"^]  E^  + 


-U   -,   /T  L  ,  n  aH-I    „  n-1 


.  [1  .  ^u  C  e""-l  .(l^AuO—  -e— 3  3^ 


1+  uC-e^^ 
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Here   E   is  the  maxlimjm  absolute  "initial"  error  which  by 
o 

(5«3)  and  (3.l5)  vanishes  indentically.   However,  we  retain 
this  term  to  aide  in  the  stability  proof. 

We  now  take  the  limit  in  (S.li;)  as  A\l   — >  0  while 
£iU  satisfies  (^•2)»  Then  the  lethargy  mesh  becomes  and 
remains  so  fine  that  condition  (i|.5)  is  satisfied  and  a  >  0 
(as  discussed  in  Section  I|.)»   Thus  by  the  definitions  in 
(5»11),  for  sufficiently  small   Au,   B   is  independent  of 
the  mesh  spacing  and  the  limit  yields,  since  ^u   =  u  /n, 

(C-l)u 
(5.15)    lim  E  <  [Ce      ^  lim  { An  E   )  + 

AU->0  ^^u-s»■0 

p     (C-l)u„  , 
+  [1  +  C   (g      n-1  _  ^3  g   lim  T  . 

^u->0 
By  (5.2),  (5.3)  and  (5»8)  the  above  implies  convergence 
provided  lim  I^=  0.   However,  this  is  just  the  condition 
required  for  consistency,  as  may  be  seen  from  (5«0),  (5»1) 
and  (5«5)»  and  it  is  satisfied  since  the  mesh  was  chosen  to 
satisfy  (5»2)»  The  convergence  proof  is  thus  complete. 

The  stability  proof  is  almost  identical  with  the 
convergence  proof  and  so  will  only  be  summarized  here.   Let 
\J»^(x)  be  the  solution  of  the  system  (3»ll)  and  (3«li4-)  when 

J 

^  .(x)   is  replaced  by   ,^  .(x).   Then  the  approximating 
system  is  said  to  be  stable  if 


-  26  - 


(5.16)  ^l/^(x)  -^  ^^(x)   as   ^j(x)  ->  ^^(x)  , 
Introducing 

(5.17)  P„  s  max   Uj(x)  -  jli^{x)\       , 

^   -a<x<a  ^       ^ 

-J<k<J 
we  find  as  above  that   P   is  bounded  by  the  expression 
on  the  right  hand  side  of  (5»ll4-)  where  now 


(5.18)  '^   =   max 
x,n,k 


Jk(^)  -=ik(^^ 


E  =  B  max 
°     x,k 


lS(x)  -i°(x) 


On  any  fixed  mesh,  the  coefficients  of  E   and   t:  are 
bounded  and  hence  (5«l6)  follows  immediately.   As  the  mesh 
is  refined  these  coefficients  remain  bounded,  as  shown  in 
deriving  (5.l5),  and  hence  the  stability  is  demonstrated  on 
all  sequences  of  meshes  for  which  the  condition  (l4.s5)  becomes 
and  remains  valid.   This  includes  all  meshes  satisfying  {S,2) 
for  which  the  approximating  equations  have  been  shown  to  be 
consistent  and  convergent.   It  can  be  shown^,  however,  that 
there  are  sequences  of  meshes  on  which  stability  holds  but 
not  consistency  (e.g.  a  sequence  on  which   A(j,  =  [^u]  ), 


^  By  more  careful  arguments  an  equality  can  be  obtained 
in  place  of  (5.1 )«  Then  a  condition  similar  to  {S»2) 
is  obtained  as  a  necessary  condition  for  consistency. 
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PART  II 

Monoenergetlc  Isotropic  Scattering 

6«  Formulation  of  Isotropic  Scattering  Problems, 

For  monoenergetlc  (one  velocity)  neutrons  which  scatter 
isotroplcally  in  a  homogeneous,  plane,  medium  the  steady 
state  transport  equation  (2o0)  becomes 


(6,0a) 


)  ^^  ix  "^  ^T  ( "i^^'^^   "  2fi^  =  J(x,ti)  , 


where   the   scattering  terra   is  now 

1 
(6.0b)  S[5]   =  o^  I    vT  I(x,ti)  dti     , 

Here   o"™   and   c   are  constants  and  c   is  the  average 
number  of  neutrons  produced  for  each  neutron  which  suffers 
a  collision.   Thus  if  no  fissionable  material  is  present 

(6.1a)  0  <  c  =  -r  <  1   , 

o-^  -    ' 

where     o"      is  the  macroscopic   scattering  cross-section. 
If  fission  may   occur  then 

(6,1b)  0  <  c   =     ^^    ^        , 

where  <T^     is  the  macroscopic  fission  cross-section  and 
V   is  the  average  number  of  neutrons  emitted  per  fission. 
In  this  case  the  medium  may  be  chosen  such  that   c  >  1, 
As  in  Section  2  boundary  conditions  are  obtained  by 
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specifying  the  flux  of  neutrons  incident  on  the  slab, 
-a  <  X  <  a*   If  there  are  known  sources  present  in  the 
vacuiAin  external  to  the  slab  the  boundary  conditions  are 
of  the  form 

(6.2)     i(-a,ti)  =  i(M-),   0  <  ti  <  1;  5(a,ix)  =  rC^i),  -l<[i.<0 

Here   X(p.)   and   r(|j,)   are  given  functions  which  determine 
the  incident  angular  flux  distribution  from  the  left  and 
right  respectively. 

The  above  formulation  includes  a  variety  of  special 
cases  which  are  of  interest  in  reactor  theory.   We  shall 
formulate  below  three  of  the  more  important  problems.  The 
approximate  solution  of  these  problems  is  considered  in 
Sections  8-10. 
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6Ao   Reactor  Crltlcalltj. 

Let  the  slab  be  composed  of  fissionable  material  for 
which  c  >  1.  Let  there  be  no  inhomogeneous  sources  of 
neutrons  present}  that  is 

(6.3)        jiif(x,vt)  =  0  ,     Ji{\L)   =  r(ti)  =  0  . 

The  criticality  problem  is  to  find  the  relationship  between 
c,  CTrn     and  a  for  which  (6o0,,2,.3)  has  a  finite  non-zero 
solution  which  is  of  one  sign,  say  positive.   This  is  an 
eigenvalue  problem  in  which  any  two  of  the  quantities  (c, 
Ow»  a)  may  be  given  and  the  third  is  then  the  eigenvalue 
to  be  detennined  such  that  the  corresponding  eigenfunction 
|(x,ii)  >  0. 

Physically  the  problem  is  to  adjust  the  material  constants 
c   and  <SL     with  the  thickness  of  the  slab,   2a,   so  that 
the  number  of  neutrons  escaping  from  the  slab  is  just  equal 
to  the  excess  of  those  produced  over  those  absorbed  in  the 
slab. 


Physically  it  is  obvious  that  more  than  one  neutron 
must  be  produced,  on  the  average,  per  collision  to 
sustain  a  chain  reaction. 
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6Bo      Capture  Fractlorio 

Let      c   <   1,      (SZ,      and     a     be  given  and    assiune  no 
distributed  scl-jrces   present   in  the   slab}    that   is 


(6.1;)  y^{^,\i)   =  0    . 

Let  the  incident  flujc  be  specified,  as  in  {6<,2)o   Then 
the  number  of  neutrons  entering  the  slab  per  unit  area 

i^  o  1 


(6.5a) 


-1  o 


1 

=  S  V'   ['^(t^)  +  r(-^i)]  d^  , 
o 

The  number  of  neutrons  captxared  in  a  cylinder  of  unit 

cross-sectional  area  through  the  slab  and  normal  to  its 

sides  is 

a     1 

(6»5b)         N  =  cr™(l-c)  j  dx     v)  d^i  i(x,M.)  , 

°    ^      -a    -1 

where  J(x,(j,)   is  the  solution  of  ( 6.0,,2,  oi;) .   Then  the 

fraction  of  neutrons  captured  in  the  slab  is 

(6.6)  P  =  /  ° 

e 

The  variation  of  P  with  c,  cJl,  a   and  the  incident 

flux  distribution  is  desiredo 
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6c «,   Escape  Probability. 

Let  no  neutrons  enter  the  slab  from  external  sources  J 
then 

(6o7)  i(ti)  5  r(ti)  =  0  • 

The  number  of  neutrons  Injected  by  the  prescribed  distributed 

source,   ^^'(xjH,),   Into  a  cylinder  of  unit  cross-sectional  area 

through  the  slab  and  normal  to  Its  sides  is 

a     1 

(6.8a)  N  =  f  dx  J  dy,^(x,ti)  . 

°   -a    -1 

The  number  of  neutrons  escaping  from  a  unit  area  on  the  sides 

of  the  slab,   x  =  +  a,   is 


ix  J^Q-sM-)  dti,  f 


K  =  -  i    V'   i(-a,M.)  dti  + 

^    -1 

(6.8b)  ^ 

=  J  |x  [i(a,ix)  +  i("-a,-^L)]  d^L  , 
o 

Here  J(x,ii)   is  the  solution  of  (6o0,o2,o7)o   The  fraction 

of  source  neutrons  which  escape  from  the  slab,  or  the  escape 

probability,  is  then 

^x 
(6.9)  ^  =  iT  • 

o 

The  variation  of  P  with  the  given  parameters   c  <  1,  cC,   and 

a   is  desired.   In  escape  probability  problems  the  source  is 

usually  taken  to  be  Independent  of  position  and  frequently 

isotropic. 
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7»   Approximating  Systeina«» 

As  in  Section  3  of  Part  I  we  shall  replace  the  continuous 
variable   p,  by  a  set  of   2 J  mesh  points,   p,..   This  procedure, 
when  applied  to  the  present  transport  problems,  is  just  the 
¥ick-Chandrasekhar  method  [10,11],   Most  applications  of  this 
method  which  appear  in  the  literature  are  to  infinite  or  semi- 
infinite  regions.   However,  we  consider  here  only  slabs  of 
finite  thickness.   In  addition  the  present  treatment  of  the 
method  is  somewhat  more  general  than  that  usually  given  and 
the  resulting  algebraic  problems  are  reduced  to  systems  or 
matrices  of  order   J  rather  than  of  order  2J<,   In  fact  the 
solution  of  a  large  class  of  problems  is  reduced  to  the 
inversion  of  at  most  two  or  three  matrices  of  order   J. 

The   2J  mesh  points  are  required  to  satisfy 

(7.0)  0  <  \i^  <  \x^  <    o,e   <  ^ij  <   1    i        tx_.  =  -|i,  ,    1  1  J  1  J- 

At  each  point  of  this  mesh  the  flux  and  scattering  term  are 
approximated  by 

(7.1)  i(x,p,,)  ^    ct.(x)  ,     S[2(x,ii^J]  ^  S[(t^(x)]  o 
Setting  ^x  =  n, .   in  (6,0a)   and  using  these  approximations 

J 


(7.2a)  ,  ^     ^ 


J 


Here   ;?3.(x)  s  )i  i^sii^)      and  tne  form  of  the  approximate 
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scattering  term  is  obtained  by  using  a  niimerical  quadrature 
formula  in  (6«0b)o  By  the  symmetry  of  the  mesh  such  an 
approximation  can  be  written  as 

(7.2b)      S[4j^(x)]  s  6-^  I  21:  OJi,  [<ti,(x)  +  (t.i^(x)]  , 

The  CO,       are  the  weights  for  a  numerical  integration  over 
the  interval   0  <  |i  <  1  using  the  points   p,,>0  of(7oO)» 
If  we  require  that  the  formula  be  exact  for  polynomials  in 
\i     of  degree  <  N  these  weights  must  satisfy 

J       n    1 
(7o3)      "^k  ^  °   5  ^  ^k  he  "^  H+1  »    n  =  0,1,2, .oe,N, 

By  a  proper  choice  of  the  \x .      and   CO.  we  may  have 

J  J 

N  =  2J-lj     but  no  greater   [8,10] «      For  the  present  analysis 
we   require   only     N  =  lo 

The  boundary  conditions    (6o2)    are    replaced  by 

il.h)  4j(-a)  =    -^j    p      1  <   J  <   JS     4j(a)   =  r^    ,   -J5J5-I5 

where   /=  =  J?(ixJ  and  r,  =  r(iio).   Equations  (7<.2)  and 
(7«i4.)  are  the  approximating  systems  v/hose  solutions  are 
sought • 

For  the  determination  and  analysis  of  these  solutions 
it  is  convenient  to  formulate  the  approximating  system  in 
matrix  notation*   For  this  purpose  we  introduce?   the 
J-dimensional  column  vectors 


3k. 


(7.5a)   ^^ix)   = 


^K.i^)\ 


Yt+j(x) 


,i,(x)   H 


/Xl(-)^     f^ 


f  ^-l' 


.    L   = 


\J+jiy^)  j 


.1 


,    R   = 


•       » 


J/ 


the   J-order  square  matrices 


0 


(7.5b)     m  = 


^2 


\0  *    ^^j/ 


»  JT  = 


^1^2    • • • 

"j 

0              « 

• 

0                « 

• 

p  =  m"      -   q. 


^q  =  I  m'^  Jli 


the  2J-diinensional    column  vectors 

4^.(x) 


(7.5c)  <Kx)  = 


J(^)   s 


(fjx)  / 


>i4(x) 


and  the   2 J-order   square  matrices S 
(7.5d) 


m    0 
M  =  [         ,  , 
0   -m 


A  = 


P   -q 


\ 


q    -p  / 

Using  these  definitions  in  (7.2)  we  obtain  the  matrix 
equation 


(7.6) 


^  (t(x)  +  CJJ  A  (|.(x)  =  K'^JU) 


The  boundary  conditions  (7.U-)  become 


(7.7) 


<^^(-a)  =  L   , 


(^  (a)  =  R 


By  using   the   same   niimerical    integration  formula   in 
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(6,5)  and  (6.8)  as  is  used  to  obtain  (7«2b)  we  may  define 
approximations  to  the  capture  fraction,   f^  a  P,   and 
escape  probability  p   «5  p,   which  are  "consistent"  with 
the  approximation  of  the  fliix.  To  employ  the  notation  of 
(7»5)  In  these  quantities  we  introduce  the  J-order  row 
vector,  which  is  a  row  of  _f^  , 

(7.8)  (J  =   (  (^1,^2*  •••»^j)  • 

The  quantities  (6,5)  and  (6»8)  are  then  approximated  by 

a)  N  J5^  n^  =  U)  m[L+R]  , 

©     e     ^^ 

a 

b)  N  ^  n^  =  <5;;(l-c)  (U  S     U+(x)+  <t  (x)3  dx  , 

>  -a 

(7»9 )  ^ 

c)  Nq  »  n^  =  (^     [i^.(x)+J'_(x)]  dx  , 

-a 

d)  N^  «  n  =  U)  ra[<t.(a)  +  <i>_(-a)]  « 

The  approximate  capture  fraction  and  escape  probability  are 
given  by 

(7.10)       fj^^  »        Pj^H^  • 

e  o 
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8>  Analysis  of  Formal  Solutions, 

The  system  (7<»6)  is  a  system  of  first  order  ordinary 
differential  equations  with  constant  coefficients.   Solutions 
of  such  a  system  can  be  obtained,  in  principle,  by  the 
application  of  well  known  formal  procedures  [12],   We  shall 
employ  the  methods  of  matrix  calculus  [13] •  Thus  w© 
introduce 

(8.0)  e'-^y  E  B(y)  =  ' 

B2i(y)  B22(y) 

where      B(y)      is    a  2J-order  matrix  and   the     B.  .(y)      are 
J-order  matrices.      The   general   solution  of    (7 96)   can  then 
be  written  as 

J 

(8.1)  <Kx)   =  B(y)5   +  ;^      r     B(y-y')    M"^;j'(y'  )    dy'     . 

^T      o 

Here  we  have  introduced  the  dimensionless  variable 

(6.2)  y  =  a;^x,    -^  <  J  <  ^   ,       bscr^a, 

which  measures   distance   in  units   of   the    total  mean-free-path, 
l/cr^    ,      of  the    slab  material.      The   2J-dimensional   column 
vector 


(8.3)  4  E 


;+i 


is  the  flux  at  the  center  of  the  slab  and  can  be  determined 
by  using  (8.1)  in  the  boundary  conditions  (7,7),   However, 
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before  treating  particular  problems  we  shall  derive  some 
properties  of  the  matrices  A  and  B(y)  which  are  required 
for  the  applicationSo 

Let  the  powers  of  A  be  denoted  by 


(8oU) 


in 


An)        An) 
^11    *12 


.(n)    (n) 
^21    *22 


where  the  a) ,   are  J-order  matrices«  Then  we  have 
Lemma  I»        A^J^  =  (-1)^  A^^^  ,   A^^^  =  (-1)^  A^J^  . 
Proof.   Since  A°  =  I  ^   the  2J-order  identity,  the  result 
is  true  for  n  =  0,   (Also  by  the  definition  (7<,5d)  of  A 
it  is  true  for  n  =  lo)  For  an  induction  we  assume  the 
lemma  true  up  to  some  arbitrary  n«  Then 

(8.5)   a"*i  =  a  a"  = 

[qA^f    -  pA<j'l    UA<|>    -  pa(|)) 

and  by  means   of   the   inductive  hypothesis 

4r'   -  P^t'   -   '^4f   =   <-!'"  tpA<;'    -  qA(j']    -   (-D-lA^fH 

Similarly  it  is  seen  that  A  j^"^-""  ^  =  (-l)^"^^  ^zT^^      ^^  ^^^ 
proof  is  completed. 

Prom  the  power  series  which  defines  the  exponential  of  a 
matrix  and  definitions  (8oO)  and  (8oU)  we  have 
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CD  /   \n   /  \ 

(8.6)  B..(y)  =YZ^^=^^^[V    • 

^J      n=0   nl    ^'^ 

Using  Lemma  I  and  these  expressions  we  obtain 

Lemma  II;        ^11  ^^^  ^  ^22^"'^^  '   ^12^^^  "^  ^^^^(-y)  . 
By  means  of  this  lemma  the  evaluation  of  B(y)   is  reduced 
to  the  evaluation  of  only  tvo  J-order  matrices. 

Another  useful  result  is  contained  in 
Lemma  III;       For  all  n  >  1  , 

4l^  ±  ^2^   =  [P  +  (-D'^'^q]  [P  +  {-1)'"'\]    ...[p  +  q]. 

Proof.  Prom  the  definitions  (7«5d)  and  (8.1|.)  the  lemma  is 
true  for  n  =  la  For  an  induction  we  assume  it  to  be  true 
up  to  n.  Then  from  (8.5) 

=  [p  +  (-1)''  q]  U[l^   ±^2^]       , 

where  we  have  used  Lemma  I.   The  lemma  now  follows  by  an 
application  of  the  inductive  hypothesis. 

Using  Lemma  III  and  the  definitions  (8,6)  we  may  write 

00   ^2n 

(8.7)  B,,(y)+B,p(y)  =  JZ  -^ [  (p+q)  (p  +  q)]""  - 

■'"'-   ~  -^^      n=0  (2n)l 

CD    2n+l 

-  (P  +  q)  H  -^ r  t(p±q)(p  +  q)]""  . 

n=0  (2n+l)L 
As  will  be  shown  in  the  next  section  it  is  the  above  two 
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j-order  matrices  which  determine  the  solution  of  all  problems 
with  homogeneous  distributed  sourceso   To  obtain  a  simpler 
representation  of  these  matrices  we  must  examine  the  matrix 

(8c8)  K^  =  (p+q)(p-q)  =  m"^  (I  -  c  JX_)    . 

Here  we  have  used  the  definition  (7»5b)  and  the  justification 
for  introducing  the  square  of  a  matrix,   Kj   is  as  follows? 
AssTjme  K   to  have  simple  elementary  divisors   (or  what  is 
sufficient,  that  its  eigenvalues  are  distinctg  see  Theorem  I ) « 
Then  there  exists  a  diagonalizing  similarity  transformation  [13] 

(8.9)  P"^  K^  P  =  A  =  i\   5^j)  o 
If  we  define 

(8.10)  k^  5V^  p  A^       =  ^\   ^ij)  > 
then  the  matrix  K  may  be  chosen  as 

(8.11)  K  =  P  /^^^   P""*-  , 

This  is  a  standard  procedure  for  defining  the  "square-root* 

of  a  matrix.   The  columns  of  the  diagonalizing  matrix   P 

2 
are  proportional  to  the  eigenvectors  of  K  . 

An  analysis  of  the  eigenvectors  and  eigenvalues  of 

2 

K   is  contained  in  /  u  \ 

Theorem  I?   (a)  The  eigenvector  u  =  1  »   I  belonging  to  an 
eigenvalue  \     of  K   has  components 


1^0 


( o » 12 ;  u  » — — 5  « 


2 
(b )   The  eigenvalues   X .   of  K   are  the   J 

J 


roots  of 


(8,13)        Pj(X)  sZI— ^  =^  • 

•^     k=l  1-Xti^   ^ 

(c)  The  eigenvalues   X.   are  real  and  distinct 
and  lie  In  the  Intervals  I 

^^J  ^^J"1         ^^j+1 

2    ••"    -^1    2 

(d)  All  eigenvalues  are  monotone  decreasing 
functions  of  c   and  the  smallest  one  has  the  values; 


(8.1S)  .J 


C  =  00 


I  "00  , 

Proof?   An  eigenvector  u  and  its  corresponding  eigenvalue 
X  must  satisfy 

(8.16)  K^  u  =  X  u  , 

From  the  definitions  (6,8)  and  (7»5b)  we  find  that  the  j-th 
component  of  (8 ,16)  Is 

-2         ^ 

-  ill  » 


and  thus 


CJ 


k  "k 


^J  = 


1-X  ixj 


Since  the  numerator  Is  independent  of   j   part  (a)  of  the 
theorem  follows.   Using  the  components  (8ol2)  of  an  eigen- 
vector in  (8.16)  we  find,  as  above,  that  the  corresponding 
eigenvalue  must  satisfy  (8813  )o   This  concludes  the  proof 
of  part  (b). 

To  examine  the  roots  of  (8el3)  we  observe,  by  (7»3), 
that   ^,  >  0  and  hence  the  k-th  term  in  the  sum  is  a 

piecewise-monotone  increasing  function  of  ?^  with  an 

-2 
infinite  discontinuity  at  X  =  [ly.   f,     The  function  Pt(^) 

is  thus  composed  of  J+1  branches  as  indicated  in  Figure  1. 

J 
Since  J     ' U),    =  1,  by  (7.3),  we  have  Pt.(0)  =  1.  The  roots 


Figure  1.   Schematic  graph  of  Fj(?^)  with   J  =  5. 
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of    (8,13)   are   the   intersections    of   this   curve  with  the 
horizontal   line  with  ordinate     l/c »      Parts    (c)   and    (d)    of 
the   theorem  are   iimnediately   obvious   from  Figure   1  and    the 
proof   is   complete. 

Prom  the  definitions    (8,8)    and    (7o5t>)   we   have 

(8.16)  (p+q)    =  m"^    ,         (p-q)    =  m  K^    ,       (p-q)(p+q)    =  ra  K^  m'^. 

Using  these  expressions  and  (bo8)  in  (8.7)  we  obtain 

a)      B^^(y)+B^(y)    =   [cosh    (Ky)    -    (mK)    sinh(Ky)3    , 


(8,1?) 


b)      B^^(y)-B^2(y^    =  xn[cosh    (Ky)    -    (Km)'^  sinh(Ey )  ]m"^. 


The  definitions    of   the  hyperbolic  functions    of   a  matrix  are 
of  course   the   series    expansions.      In    (8.17b)   we   have  used 
K"       which  exists    only  if     c   7^  1.      The   explicit  dependence 
on     y      of    the  above  matrices   may   be   determined  by  using 
(8.10-„11)    in    (8,17)    to   obtain 

a)  ^ii(j)-^^i2^j)   =  ^^^^ij   cosh  k^.y) 

-  (P"^mP)(6^^.kj   sinh  k^y)]   P"^   , 
(8.18) 

b)  B^-|^(y)=B^2^3r)    =    ^^^^    [  (^ij   cosh  k^.y)    - 

-  (P"^  m°^   ^^^^i1   ~¥T  ^^^^  k^y)](mP)"^. 

Another  matrix  of   interest   is     A"    <,      To   obtain  this 
matrix  we  first   observe   that  by    (7.5^)    and    (7«.3)   with     n  =   0: 


i+3  - 


2 

JX     =   Sl  ' 

Then  from  the  definition    (7.5d)    of     A     we  can  easily  verify- 
that 

r       3\  Vr=m-s 

c 


(8.19)  A.'^  =\  I      I 


-3     -r 


^   =  2{c-l)  -Hldi  i 


provided      c   7^  1«      The   singular   case,      c   =  1,      will  not   be 
treated   here   in  any  detail  as   the   present  results    hold  for 
c     arbitrarily  close   to   one.      However  this  exceptional   case 
offers   no  formal  difficulties  and   Indeed  is    the   first  problem 
studied,  for  semi-infinite   slabs,    by  Chandrasekhar    [10]« 
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9o      General  Solution  for   Homogeneous   Sotirces. 

The  general   problem  formulated    in    (7 •6-. 7)    can  be 
solved   explicitly  if   the  source,      j6)(x),    is    independent 
of  position,,      We  consider   here  such  problems    in  which 

(9oO)  W^(x)   =  siJ  (=  a  constant   vector)    • 

The  general  solution    (8»1)   now  becomes 

y 

(9<.l)  <t>(x)    =  B(y)    «   +  ^     J'     B(y-y')    dy«    M'^     . 

°T      o 

From  the  definition    (8.0)    of     B(y)      we   have 

J 
(9c2)  J'B(»y')   dyt    =  A"^  B(-y)    =  B(-y)   A"^    , 

since   B(y)   is  a  power  series  in  A.   Using  (9«2)  in  (9ol) 
we  may  write  the  solution  as 

<t>(x)  =  B(y)  e  +  ^  [I-B(y)]  A'^  n"^  J     , 

T 

(9.3)  ^ 

=  B(y)  (^-  ^)  +  ^^  , 


where  we   have   introduced   the   2J-dimensional   vector 

/''^^  1      -1-15 

The  boundary  conditions  (7«7),  from  which  C  is  to  be 
determined,  become  by  (9 •3) 

<t)^(.a)  =  B^^(.b)(C+»)7+)  +  B^2("^)(^.-^/-)  +^+  =  L  , 

(9.^) 

<|)^(a)   =  B2-L(b)(C+-^4.)  +  B22(b)(C.-y„)  +  >^.  =  R  • 
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Using  Lemma  II  the  solution  of  this  system  can  be  written 
as 

(9.6)  4+  =  »!+  +  I  [B22(b)+B2-L(b)]-^[(L+R)-()^^+^/_)]  + 

+  1  lB^^{h)-B^^{b)r^   [(L.R)-(V|^-)^_)]  . 

The  solution  (9»3)  Is  thus  determined,  and  by  another 
application  of  Lemma  II  and  (9 •6)  it  may  be  written 
explicitly  as 


-Ix 


(9.7) 


<i+(x)  =  »]+  +  |[B22(+y)+B2i(:fy)][B22(b)+B2i(b)3 

x[(L+R)-(^^+'/_)]  + 


i|[B22(-?y)-B2i(:pr)][B22(b)-=B2-L(b)]-lx 


>t[(L-R)-(^^-^_)]  . 

As  mentioned  in  the  previous  section  we  see  that  the  solution 
depends  only  on  the  matrices  given  in  (8.1?^  or  (8,l8) 
(recalling  Lemma  II )»   The  only  matrices  which  are  not 
given  explicitly  are   P"  ,   to  be  used  in  (8,18),  and 
[Bp„(b  )+Bp,  (b)]""  .   These  are  matrices  of  order   J  and, 
for  a  particular  problem,  are  independent  of  the  spatial 
variable.   For  many  problems  of  interest,  including  isotropic 
soiarces,   L-R  =  '^+-'7_  =  0   and  hence  [B^^{h)-B^^{h)]''^     is 
not  required. 
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The  general  solution  (9<.7)  may  be  simplified  to  study 
a  variety  of  special  problemso   In  the  next  section  we 
consider  those  problems  formulated  in  Sections  6Aj,  B  and  Co 
The  study  of  the  reflection  and  transmission  matrices  of 
slabs  J  which  are  of  more  interest  in  astrophysical  problems 
[10],  offers  no  difficulty  but  will  be  omitted  here* 
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10.   Special  Problems, 

lOAc   Determination  of  Criticality. 

An  approximate  solution  of  the  criticality  problem 
fonnulated  in  Section  6A  is  obtained  by  using  (6,3)  in 
(7«>5a).   Then  we  have 

(10.1)         ;>^(x)  SO,         L  =  R  =  0  , 

and  the  formal  solution  (9«>1)  becomes 


(10.2) 


4(x)  =  B(y)  4 


The  boundary  conditions  (7»7)  are,  by  (10.1)  and  Lemma  II, 


(10.3) 


B^^ih)   4^  +  B23^(b)  4.  =  0  , 
B2i(b)  ?^  +  B22(b)  C„  =  0  . 


These  equations  are  homogeneous,  in  contrast  to  those  of 
(9.5)»  and  have  a  solution  only  if  the  coefficient  matrix 
is  singular,  that  is 


(lO.ij.) 


B^^{h)        B2i(b) 


B2]^(b)   B22(b) 


=  0 


The  solution  of  Section  9  is  thus  not  applicable  and  we  are 
faced  with  the  problem  of  find: 
for  which  (10. U)  is  satisfied. 


faced  with  the  problem  of  finding  values  of  c   and  b  =  cC  a 


-  kQ  - 


However,  the  above  condition  may  be  simplified  by 
observing  that  the  solution  of  interest  must  have  the 
symmetry  property 

(10.5)  <t+(x)  =  4„(-x)  . 

This  is  physically  obvious  for  the  solution  of  the  exact 
problem  from  the  symmetry  of  that  problem^   More  precisely 
if  5(^»i^)   ^^  ^^  eigenfunction  then  so  is  J(-x,-^l),   as 
may  be  verified  by  substitution  into  (6<,0).   The  symmetry 
of  the  exact  solution  is  then  a  consequence  of  the  uniqueness 
of  the  eigenf\inctions»   A  similar  result  is  also  true  of  the 
eigenvectors  J   (i(x),   of  the  approximating  system,  from  which 
(10<,5)  follows.  From  (10«5)  we  have  in  particular 

4+(0)  =  4_(0)  , 

and  since  B(0)  =  I  we  deduce  from  (10,2)  that 

(10.6)  §+=€.« 

The  boundary  conditions  (10.3)  are  now  reduced  to 

(10.7)  f^22^^^  "^  ^21^^^^  ^+  =  °  ' 

a  homogeneous  system  of  order   J. 

Non-trivial  solutions  of  (10.7)  exist  provided  that 

(10.8)  |B22(b)  +  B2-,^(b)!  =  0  . 


k9 


Equation  (IO08)  is  called  the  criticality  equation,,  It  may 
be  determined  explicitly  by  using  Lemma  II  and  (8.l8a)<.  If 
cosh  k.b  /  0   for  all   j   the  criticality  equation  becomes 

(10.9)  1 1  +  (P'-'-mPXe^  .k  tanh  k.b)|  =0  . 

For  each  root,   (c,b),   of  (10o9)  K+     is  determined,  to 
within  a  scalar  factor,  as  a  solution  of  (10«7)o   The 
desired  solution  is  that  for  which  (|"^(x)  >  0  in  -a  <  x  <  a. 
Prom  (10<,2),  Lemma  II,  and  (10»6)  the  solution  is  found  to  be, 
formally, 

(10.10)  ^J-^)   =  <t+(x)  =  [B22(-y)+B2i(-y)3  ^+  o 
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lOBe   Calculation  of  Capture  Fraction, 

The  solution  of  the  capture  fraction  problem  formulated 
in  Section  6B  is  approximated  by  using  (S,!;)  in  (7»5a),  The 
solution  of  Section  9  now  applies  with  y3  =  0«  Then  by  {9,l\.), 
^5  0,   and  from  (9<i7)  we  obtain 

x[B22(b)+B2]^(b)]~^  (L+R)  , 

This  expression  may  be  simplified  by  using  Lemma  XI  and 
(8.17a)  to  yield 

(10.11)  (^^(x)+<i_(x)  =  [cosh  Kj][B^^{h)+B^^{h)r^    (L+R)  . 

Recalling  that  y  =  o^  x  we  obtain 
a 

(10.12)  /  [(i^(x)+(t_(x)]  dx  =  ^  K~^[sinh  Kb]  [B22(b)+B2-L(b)]"^(L+Ri 

For  simplicity  we  introduce  the  notation 

(10.13)  Q,  =  K"-"-    [sinh  Kb]    [B22(b  )+B2-L(b  )]"■'■    , 
and  use   it   in   (10.12)   and    (7.9b)   to   obtain 

(10, lU)  n„    =  2(l-c)      U)     Q(L+R)    , 

Then  from    (7«9a)   and    (7. 10)   the   capture   fraction  is   fo\and 
to  be 


(10.15)  fj  =  2(l-c) 


W   i^(L+R) 
Cj   m(L+R) 
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lOCo   Calculation  of  Escape  Probabllltyo 

The  escape  probability  problem  formulated  in  Section  6C 
is  approximated  by  using  (6o7)  In  (7»5a)  to  get   L  =  R  =  Oo 
Taking  the  distributed  source  to  be  independent  of  position 
the  solution  of  Section  9  applies  and  from  it  we  obtain 

4^.(a)+<|._(-a)  =  {^^+>lJ~[B^^{-h)+B^^{-h)][B^^{h)+B^^{h)r\yj^+^J, 

Using  Lemma  II  and  (8ol7a)  this  expression  may  be  written  as 

4^.(a)-»4__(-a)  =  2inK[sinh  Kb][B^^{\>)+B^^{h)r^    (7++7->  * 

(10,16)       =  2m  K^  Q  i^+*^J    » 

where  we  have  used  the  definition  (10«13)»   Prom  (9oi4.)j  (8o8) 
and  (8«19)  we  obtain 

(10.17) 

Here  we  have  used  the  easily  verified  result? 

(I  -  c  .TL)"-"-  =  (I  -  -^  _rL)»  Using  (10,17)  and  (10,16) 

in  (7,9d)  we  have 

(10,18)        ^  =  ^  '^>  (^^  K^)  <^(^^  K^)""^  (?^+  +^_)« 

Since  the  source  was  assumed  independent  of  position 
(7.9c)  gives 
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(10.19)  n^   =  2a  oJ    (/^+  +  4  J    » 

and  the  escape   probability  Is,   from    (TolO) 

6J   (mV)^(inV)"^(i^+;^_) 


(10,20)        Pt  =  —r     — n        /j — 


This  fonm  of  Pj  has  been  derived  to  show  some  formal 
relation  with  the  expression  (10.15)  for  f -.  A  better 
representation  for  computations  is  derived  in  the  next 
section}  equation  (10.26). 
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IQDo   A  Relation  between  Escape  Probability  and  Capture  Fractlono 

An  Interesting  relation  between  the  escape  probability 
and  capture  fraction  can  be  obtained  when  the  angular 
dependence  of  the  respective  source  terms  are  related  by 

(10.21)  (xf+  +  xf _)  =  constoX(L+R)  s  S  « 

In  terras  of  the  continuous  source  functions  this  implies  that 

(10.22)  [J{\i)+j{-\i)]    =  const o U( IX )+r(-tjL)]  =  s(ijl),   0  <  ^l  <  1. 

Here   s(tji)   Is  a  function  Introduced  for  convenience  of 
notation  and  the  J-dlmensional  column  vector  S  has 
components   S«  =  s([i.,). 

To  obtain  the  desired  relation  we  first  simplify  (10.20) 
by  noting,  from  (8.8),  (7.5b),  (7.8)  and  (7.3)  with  n  =  0, 
that 

6J  (ra%^)  =  U)  (I  -  c  _n.)  , 

(10,23) 

=  (1"C)  CO 

Prom   (10.17)  we   have,   using  the   notation   (10.21), 
(10.21+)  (mV)"^^++i.)   =  S  +  3^  [CaJ'S]    1^. 

We  note   that      [  Cjl)  »  S]      is    a   scalar  and  we  have   Introduced 
the   J-dlmenslonal  coliomn  vector 

h 


(10.25)  1   s 


u 


5U 


Using    (10.21),    (10.23)    and    (10.2lj.)    in    (10.20)  we   obtain 

(10.26)  Pt(S)=|;|     ^Ts""^^      UJ      Q  1    , 

Here  we  have  included  the  argixment  S  in  p^  to  explicitly 
indicate  the  dependence  on  the  angular  distribution  of  the 
source. 

With  similar  notation,  from  (10.21),  the  capture 
fraction  (10.15)  can  be  written  as 

(10.27)  "^  ^ 

=  2(1"C) 


^J 


(S)   ^2 


We  have  introduced  here  the  source-weighted  average  cosine 
(10,28)(a)     ilj(S)  =  ^~-g-  . 

Prom  the  definitions  of  (a)  ,   m  and   S   it  is  seen  that 

this  is  a  numeric  ail  approximation  of 

1 
_      J'p,s(p,)d(j,    _ 
(10.28)(b)     p,(s)  =  J— ^   tij(S)  . 

J^s((i)dp, 

Furthermore  by  (7.3)  we  see  that  M^j^S)  =  |I(s)  if  s(ti)  is 
a  polynomial  in  p,  of  degree  <  N-1.  Of  course  as  J  —*»  oo 
in  an  appropriate  manner  ii,j(S)  ■— s»  ^^(s)  for  all  integrable 
functions   s(^). 

An  Isotropic  source  is  represented  by  setting 
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S     =     1      ,  (laeo  3(li)     =     1)      c 

Then  from   (10o27)  for   an  isotropic    source,    and  recallirig   by 
(7.3)  that      Ca^»  1_.  =  1,     we  have 

f   (1)   =  2Slz2A      uj     Q  1   , 
Using  this   result    and    (10«27)   in   (10o26)  we  have 


lij(S) 


1^.t(1) 


frd.) 


(10,29)     pj(S)    =^fj(S)   +3^     ^^^     .J..,, 

As  this  relation  holds  for  all  meshes  (7oO)  (ioe.  for  all   J) 
we  may  assvime,  following  Chandrasekharj  that  it  is  an  exact 
relation  between  the  actual  escape  probability  and  capture 
fractions «   That  is,  in  an  obvious  notation,  we  have 


(10.30) 


The  relation  (10o30)  may  be  verified  in  two  special 
cases o   First,  if  the  incident  and  distributed  sources  are 
isotropic  then   s(n)  =  1   and  from  (10o28)  jl(l)  =  l/2  , 
and  (10.30)  becomes 


P(l)  = 


F(l) 


(10.31) 


2(l-c)(3^2a 

F(l) 
2(o^-ag)2a 


Here  we  have  used  (6»la)  and  the  result,  (10.31),  is  well 
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known  [ll4.]»   In  the  second  case  consider  purely  absorbing 

media  in  which  c  =  0   and  az,   =  cr  .      the  absorption  cross 

section.   Then  (10o30)  becomes  simply 

iKs) 
(10,32)  P(s)  =  ^r^  P(s)  o 

a 

This  result  Is  generally  known  for  arbitrary  convex  bodies 

provided  the  sources  are  isotropic  [151 o   In  the  form  (10o32) 

we  have  a  generalization  to  arbitrary  sources  but  which  is 

valid  only  for  slabs© 
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